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ABSTRACT 

The magnetorotational instability (MRI) can be a powerful mechanism amplifying the mag¬ 
netic field in core collapse supemovae. Whether initially weak magnetic fields can be ampli¬ 
fied by this instability to dynamically relevant strengths is still a matter of debate. One of the 
main uncertainties concerns the process that terminates the growth of the instability. Parasitic 
instabilities of both Kelvin-Helmholtz and tearing-mode type have been suggested to play 
a crucial role in this process, disrupting MRI channel fiows and quenching magnetic field 
amplification. We perform two-dimensional and three-dimensional sheering-disc simulations 
of a differentially rotating proto-neutron star layer in non-ideal magnetohydrodynamics with 
unprecedented high numerical accuracy, finding that Kelvin-Helmholtz parasitic modes dom¬ 
inate tearing modes in the regime of large hydrodynamic and magnetic Reynolds numbers, 
as encountered close to the surface of proto-neutron stars. They also determine the maximum 
magnetic field stress achievable during the exponential growth of the MRI. Our results are 
consistent with the theory of parasitic instabilities based on a local stability analysis. To simu¬ 
late the Kelvin-Helmholtz instabilities properly a very high numerical resolution is necessary. 
Using 9th order spatial reconstruction schemes, we find that at least 8 grid zones per MRI 
channel are necessary to simulate the growth phase of the MRI and reach an accuracy of 
~ 10% in the growth rate, while more than ~ 60 zones per channel are required to achieve 
convergent results for the value of the magnetic stress at MRI termination. 

Key words: accretion, accretion discs - MHD - instabilities - stars: magnetic field - 
supemovae: general 


1 INTRODUCTION 

Originally discovered by |Velikhov| ( |1959| ) and |Chandrasekhar| 
(|1960|l, the magne torotational instability (MRI) was suggested by 
[Balbus & Hawley|(l991| BH91 hereafter) to be the physical mech¬ 
anism driving the redistribution of angular momentum required for 
the accretion process in Keplerian discs orbiting compact objects 
(see, e.g., [Balbus & Hawley] 1998| for a review). 

Keplerian discs have a positive radial gradient in angular mo¬ 
mentum, and therefore are linearly (Rayleigh-)stable. Purely hydro- 
dynamic perturbations are unlikely to grow to amplitudes at which 
the associated stresses can account for efficient angular-momentum 
transport. In the presence of a weak magnetic field, however, a neg¬ 
ative radial gradient in the angular velocity of the disc is MRI- 
unstable, and seed perturbations can grow exponentially on time 
scales close to the rotational period. During this phase, channel 
modes develop. Channel modes are pairs of coherent radial up- and 
downflows stacked vertically and threaded by layers of magnetic 
field of alternating radial and azimuthal polarity. In these modes, 
the magnetic tension (Maxwell stress tensor) transports angular 


momentum along the field lines from the inner parts of the disc 
outwards. 

The criterion for the MRI onset can be formulated in a rather 
simple manner, even if the thermal stratification (gradients of en¬ 
tropy or molecular weight) and non-ideal effects (viscosity, resis¬ 
tivity) are included ( |Balbus|1995||Menou et al.|2004| ). This allows 
for its application beyond Keplerian discs, in particular to proto¬ 
neutron stars (PNSs) resulting from the core-collapse of rotating 
massive stars. Simplified simulations by [Akiyama et al.j l |2003| ) 
showed that such PNSs possess regions in which the MRI can 
grow on shorter time-scales than the time between the bounce and 
the successful explosion. This finding, later confirmed in multi¬ 
dimensional models (e.g. Obergaulinger et al. 2006b jCerda-Duranj 


|et al.|2007l|Sawai et al.|2013[|Sawai & Yamada|2015| i, presents the 

possibility of generating strong magnetic fields that can tap the ro¬ 
tational energy of the core, power magnetohydrodynamics (MHD) 
turbulence ( [Masada et al.|2015] l, and become a potentially impor¬ 
tant ingredient in rapidly-rotating core-collapse supemovae (CC- 
SNe). 

How much these systems are affected by the MRI crucially de¬ 
pends on both its growth rate and on the final amplitude of the seed 
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perturbations. We can give an upper limit by assuming that the MRI 
ceases to grow once the magnetic field comes close (within a factor 
of a) to equipartition with the energy of the differential rotation. 
In CCSNe, this would correspond to dynamically important field 
strengths up to 10^^ G. Similar energetic arguments can be used 
to express the MRI-generated stresses in the framework of cr-disc 
models ( |Shakura & Sunyaev|1973| l. 

This estimate neglects possible effects quenching the MRI be¬ 
fore it reaches its maximally-allowed energy, and the effects of 
buoyancy as shown by |Guilet & Miillerl ( |2015| ), who performed 
MRI simulations in the presence of buoyancy. They showed that 
the termination amplitude (at the end of the exponential growth) 
is not necessarily correlated to the magnetic field strength in the 
turbulent state that followed. The physics of the termination of the 
MRI growth remains an active field of research with many studies 
devoted to finding the value of the a parameter for the stress tensor. 
We refer, among others, to the works of |Sano et al.|(2004]l; Sano 


|& InutsuS (|2001|l; Brandenburg] (|2005|l; Fromang & Papaloizou 

' iKnoblo ~ 


(|2007|l;|Gardiner & Stone|(|2005|l;jKnobloch & Julien|(|20()^ 

The model of parasitic instabilities by Goodman & Xu fl994| 
GX94 hereafter), further studied and developed by [Latter et al.| 
( |2009| ), |Pessah & Goodman] ( |2009| ) and |Pessah| ( |2010| ) provides 
a clear physical picture of the termination mechanism. The MRI 
channel modes are characterized by a shear layer and a current 
sheet in the vertical profiles of velocity and magnetic field, respec¬ 
tively. Hence, the (laminar) channel flows can be unstable against 
secondary (or parasitic) instabilities of Kelvin-Helmholtz (KH) or 
tearing-mode (TM) type [^Initially, the role of the parasites is neg¬ 
ligible, as they grow much more slowly than the MRI. However, 
since the growth rate of the secondary instabilities is proportional 
to the channel mode amplitude, it is clear that at some stage the par¬ 
asites will grow faster than the MRI channels whose growth rate 
is constant, whereas the parasites grow exponentially with time. 
Roughly at this point, the parasitic instabilities should disrupt the 
channel modes and terminate the MRI growth, marking the tran¬ 
sition to the turbulent saturation phase ( |Pessah|[2010| ). A further 
discussion of the MRI saturated state is beyond the scope of this 
paper. 

|Pessah| ( |2010| ) analytically studied the MRI termination in 
resistive-viscous MHD by solving simplified model equations for 
the evolution of the parasitic instabilities. He identified different 
parameter-space regimes where, depending on hydrodynamic and 
magnetic Reynolds numbers, either the KH instability or the TM 
is the dominant (i.e. faster developing) secondary instability. |Ober-| 
Igaulinger et al.| ( [20T^ found that the magnetic field at the surface 
of the PNS can be enhanced w.r.t. the interior regions. If this is also 
the case for rotating cores, MHD phenomena (like, e.g. the MRI) 
should be most prominent at the PNS surface. In this region, the 
Reynolds numbers are large if the surface is located above the neu- 
trinosphere, which however may not always be the case (see Fig. 10 
in |Guilet et al.|2015] l. In this paper, we only investigate the regime 
of very high Reynolds numbers in which, according to the parasitic 
model, the MRI should be terminated by the KH instability. 

The parasitic model has not been tested with direct global nu¬ 
merical simulations of the MRI in CCSNe. |Obergaulinger et aT] 
( |2009| ) found in their semi-global 2D ideal MHD simulations 
that, because of numerical resistivity, the MRI was terminated by 
TMs. In their 3D simulations the MRI was terminated by non- 
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1 20091 classifies the types of parasitic modes differently. 


axisymmetric parasitic instabilities, although a clear identification 
of their nature was not possible. 

To test the predictions of |Pessah| ( |2010| ), we performed a set of 
2D and 3D resistive-viscous MHD simulations of the MRI. Given 
initial conditions, we studied in particular the importance of non¬ 
ideal effects by varying the values of both the uniform viscosity 
and uniform resistivity, which influence the growth of KH modes 
and of TMs, respectively. In the case of the KH instability, which 
is present already in ideal hydrodynamics, viscosity and resistivity 
alter the properties of the unstable modes merely quantitatively. On 
the other hand, TM are essentially driven by resistivity, i.e. they 
do not grow in ideal MHD, while viscosity plays a minor role by 
changing the properties of the unstable modes only quantitatively. 

Since we are mainly concerned in this work with identifying 
the type of the parasitic instability that limits the growth of the MRI 
in particular models, we performed all our simulations (except for 
one control model) with a non-zero resistivity, but simulated mod¬ 
els with both zero and non-zero viscosity. We defer a more thor¬ 
ough quantitative analysis of the influence of non-ideal effects to a 
subsequent work. 

In Sec. we give the criterion for the onset of the MRI, and 
we describe the initial stage of the instability during which channel 
modes develop. Next, we discuss possible scenarios of its termina¬ 
tion, in particular, termination via the parasitic instabilities. In Sec. 
I^we describe the numerical code and the initial setup used in our 
2D and 3D simulations. We present the results of these simulations 
in Sec.|^ and summarize our findings in Sec.[^ 

2 MRI EXPONENTIAL GROWTH PHASE AND 
TERMINATION 

2.1 Physical model 

We consider flows that can be described by the equations of 
resistive-viscous (non-ideal) MHD. In the presence of an external 


gravitational potential, O, these equations read 

5,p-hV.(pv) = 0, (1) 

dt(p\) + V • (p\ <S)\ + T) = -pVO, (2) 

d,e^ + V • -h V • T -h //(b • Vb - ^Vb^)] = -pv • VO, (3) 
5^b -h V • (v 0 b - b 0 v) = y/V^b, (4) 

V • b = 0, (5) 


where v, p, 77 , and b = B/ V4^ are the fluid velocity, the density, a 
uniform resistivity, and the redefined magnetic field B, respectively. 
The total energy density, is composed of fluid and magnetic 
contributions, i.e. e^, = s + |pv^ -h |b^ with the internal energy 
density s and the gas pressure p = p(p, £,...). The stress tensor T 
is given by 

T = [P+ ib2+p(|v-^)V-v]l-b®b-pv[Vv + (Vv)^], (6) 

where I is the unit tensor, and v and ^ are the kinematic shear and 
bulk viscosity, respectively. 

2.2 Magnetorotational instability 

We study the MRI in a small portion of the rotating star at a given 
distance r from the rotation axis, embedded in a magnetic field. For 
convenience, we use cylindrical coordinates (r, 0, z), hereafter. We 
restrict our analysis to locations close to the equatorial plane (z = 0) 
and vertical perturbation wavevectors for which the MRI is known 
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to develop fastest (see, e.g. [Balbus & Hawley|1998| l. In this case, 
we can consider a differentially rotating fluid with angular velocity 
Q and linear velocity 

V = Qr0, (7) 


threaded by a uniform vertical magnetic field 

b = bQ,z ( 8 ) 

in the local perturbation analysis. Here, 0 and z are the unit vectors 
in 0 and z direction, respectively. 

With these assumptions and ignoring dissipative effects, the 
MRI instability criterion is (c.f. |Balbus|1995t 

+ rdrii^ =N^ + k^- 40.^ < 0, (9) 

P \ P ^iPj’ 

^dr(r*a\ ( 10 ) 

ri 

are the square of the Brunt-Vaisala frequency and the epicyclic fre¬ 
quency, respectively, and Y i is the adiabatic index. 

We consider an angular velocity with a radial dependence of 
the form 

n = (11) 


where 


where flo is the angular velocity at the characteristic radius ro, and 
q is the local rotational shear given by 


dlnQ 

dlnr 


( 12 ) 


The rotation profile O is quite generic for astrophysical systems, 
e.g. for ^ = 3/2, one recovers the Keplerian profile and for differen¬ 
tially rotating stars typical values of q are in the range Q < q <312. 
The corresponding epicyclic frequency (i.e. the radial oscillation 
frequency) is 


K = ^J2(2 - q)Yl, 


(13) 


which vanishes at q - 2 (Rayleigh stability criterion limit for un¬ 
magnetised rotating fluids). We also assume that the entropy and 
composition are constant within the simulated volume, and hence 
- 0. The influence of entropy gradients on the MRI has been 
studied by |Obergaulinger et ah] ( |2009| ) and will not be discussed 
here. 

If condition (|^, is fulfilled, i.e. 0 < q <2, any perturbation of 
the form (WKB ansatz) is unstable for wavenumbers 

k,<kcni=^/2q — , (14) 

^Az 

where caz = I aIp is the Alfven speed in the vertical direction. 

BH91 considered small perturbations of the velocity, Vc, and 
the magnetic field, be, in a system whose velocity and magnetic 
field are given by Eqs. ([^ and respectively, i.e. 

V = -h Vc, (15) 

b = + be. (16) 


We note that one must also introduce appropriate perturbations of 
the other thermodynamical quantities to fulfil the MHD Eqs. ^ - 
0, but we do not mention them here explicitly. 

Eor the linearized ideal MHD equations in the incompressible 


limit, BH91 found unstable solutions, which are usually referred to 
as MRI channels, 

Vc(l;^mri) = Vce’'™'(f COS 0 V + 0 sin 0 v) sin(A:^z), ( 17 ) 

= ^ce’'™'(f COS 0 /, + ^ sin ^fc) cos(A;„riZ), ( 18 ) 

where the subscript c stands for channel, r is the unit vector in r di¬ 
rection, Vc and be are the initial amplitudes, 0v and 0^ are the angles 
between the r-axis and the direction of the velocity and magnetic 
field channels, respectively. The wavenumbers correspond to 
values of fulfilling Eq. d- To simplify the notation, we define 

Vc(0 = (19) 

bM = ( 20 ) 


and for brevity, we will often drop the explicit time dependence, 
i.e. Vc = Vc(0 and be - be{t). GX94 generalized the results of 
BH91, and they showed that the MRI channels are an exact so¬ 
lution of the ideal incompressible MHD equations in the shearing 
sheet (local) approximation. This approximation consists in trans¬ 
forming the equations to a frame corotating with a fiducial fiuid 
element and linearising the rotational profile around a radius ro, 
i.e. D(r) ^ (r - ro)5rf^(^)lro- In this frame, the gravitational force 
and the centrifugal force balance each other for initially Keplerian 
accretion discs, but the Coriolis force has to be taken into account. 
In differentially rotating stars, additional pressure gradients are nec¬ 
essary to provide an equilibrium. 

In the ideal MHD limit, the MRI growth rate and the channel 
angles, 0v and 0^, are given by 


Tinri 


{2-qf + ^q\^ 


(2-q)- 2q 



( 21 ) 


and 


arctan 




■Tm 


^Tmri^O 


7T 

'^r 


( 22 ) 

(23) 


respectively. The amplitude ratio Vc/^c is a function of k^^ and 
q (cf. [Pessah & Chan||2008| PC08 hereafter). The mode with the 
wavenumber 


(24) 

\ 4 caz 

grows fastest at a rate 

Tmri — —Y1. (25) 

Note that we use capital letters in the subscripts to refer to the 
fastest-growing mode (7 mri and kum) to distinguish it from generic 
unstable modes with growth rates y^ri < 7mri- 

Eor the fastest-growing mode, the magnetic field and the ve¬ 
locity amplitudes are related by 


where cac = bel ^/p is the Alfven speed parallel to the MRI chan¬ 
nel, and the channel angles are 0v = 7t/4 and 0^ = 3;r/4. 

The properties of MRI modes change when viscosity or re¬ 
sistivity are present in the system. PC08 generalized the results 
of GX94 and showed that MRI channels (Eqs. |T7] l and |T^ ) are 
also exact solutions of the resistive-viscous incompressible MHD 
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equations in the shearing sheet approximation. They derived ex¬ 
pressions for the growth rate the amplitude ratio and the 
channel angles 0v and 0 ^ of MRI-unstable modes for arbitrary hy¬ 
drodynamic and magnetic Reynolds numbers. Following PC08, we 
define these two (dimensionless) numbers as 


R. 

Rm 


rjO. 


(27) 

(28) 


We note that in the case of a zero entropy gradient the (most gen¬ 
eral) dispersion relation of |Menou et al.| ( [2004| reduces to the one 
analysed by PC08. Since we limit our studies to isentropic models, 
we can directly apply their results. 

In resistive-viscous MHD, the channel angles, 0v and 0^, are 
given by 


0 v = arctan 
0 ^ = arctan 


+ (rmri + fcL>?)(ymri + kj^y) 

2no(rn,ri + kl^v) 

-ao{2(y„, + kl,^) + q{v - ri)kl^] 

+ (rmri + ^L'7)(rmii + 


(29) 

(30) 


when transforming Eqs. (50) and (51) of PC08 into dimensional 
units. In general, velocity channels and magnetic channels are not 
orthogonal. However, it is evident that if y = 77 (i.e. when the 
magnetic Prandtl number = R^IRq - 1), Eq. reduces to 
Eq. ( [^ , so that the shift of njl between 0v and 0^ may also hold 
in the non ideal MHD limit. 

As the expressions for the MRI growth rate ymri and the ampli¬ 
tude ratio Vc/^c are quite complex, we do not give them here (the 
interested reader will find the details in PC08). Instead, we briefiy 
discuss some of their key physical aspects. 

Both resistivity and viscosity reduce the MRI growth rate. 
They also shift both the critical wavenumber k^ni and the most un¬ 
stable wavenumber hum to lower values than those given for ideal 
MHD by Eq. ( |14| ) and ( |24| ), respectively. This behaviour can be 
readily understood, as dissipative effects are more pronounced in 
the resistive-viscous regime (as they scale oc kj). The smaller the 
Reynolds numbers, the larger the growth rate reduction and the 
larger the shifts. Moreover, the magnetic Reynolds number has a 
bigger influence than the hydrodynamic one, because the nature 
of the MRI is more “magnetic” than “hydrodynamic”. Eor any 
Reynolds number, the MRI growth rate tends to zero for 

sufficiently short wavevectors, i.e. y^ ^ 0 for ^ 0 (see Eq.[^. 

PC08 found analytic expressions for the growth rate yMRi, and 
the wavevector Z^mri of the fastest-growing mode only in some lim¬ 
iting cases, usually for very large or very small Reynolds numbers 
(e.g. for Rm ^ 1 and ^\). We note that for R^ > 10 (which 
holds for core collapse supemovae outside the neutrinosphere) the 
ideal MHD expressions given in Eqs. and determine 
the values of the wavenumber and the growth rate of the fastest- 
growing mode with relative errors < 10 % compared to the values 
obtained by the expressions of PC08. 

The high conductivity of the degenerate matter in a super¬ 
nova core implies very high magnetic Reynolds numbers. Using 
an order-of-magnitude estimate of 77 ~ 10 ""^ 

|Duncan|1993| l, Eq. ( [^ yields 

\2 

/ 79a_ 

Rr. = 10 


^cnrs ^ (Thompson & 


,/^y/ 

\10i3Gj \ 


10 ^"^ g cm ^ 


10^ s-i 


a 


(31) 

Because the magnetic Reynolds number is so large, we can safely 


neglect resistivity in the expressions for the growth rate and the 
wavelength of the MRI. 

The molecular viscosity of supernova matter (y ~ 
0.4 cm^ s "^) is only a few orders of magnitude larger than the re¬ 
sistivity, i.e. the respective hydrodynamic Reynolds number Rq ~ 
2.5 X 10^ » 1 ( [Thompson & Duncan||l99^ . In the region be¬ 
low the neutrinosphere, however, the tight coupling between neu¬ 
trinos and matter changes the situation. According to [Guilet et'H^ 
( |2015| ) the neutrino-matter interaction results in an effective viscos¬ 
ity that varies between a few 10^ cm^s“^ deep inside the PNS and 
10 ^^ cm^s“^ near the neutrinosphere, i.e. one finds 


Re = 0.1 


bpz 

1013 G 


10 i"i g cm 3 \ /103 s 1 \ / lOii^ cm^s 


P 


a 


. (32) 


The above equation implies that the expressions for the growth 
rate and for the angles 0 ^ and 0v (Eqs. |22|25| ), which are valid 
in the limit of ideal MHD, are not satisfied inside the PNS. Out¬ 
side the neutrinosphere, the interaction between matter and neutri¬ 
nos changes its character from diffusion to free streaming. In the 
free streaming regime neutrino drag damps fluctuations not like 
an effective viscosity, but like a drag force that is independent of 
the wavelength. Its impact on the MRI is therefore very different 
from the impact of a viscosity ( [Guilet et al.|2015| l. Neutrino drag 
is relevant not only outside but also significantly below the neu¬ 
trinosphere as long as the dynamics of interest is happening at a 
wavelength shorter than the neutrino mean free path, which varies 
from a metre inside the PNS to several kilometre near its surface 
(see Pig. 5 in [Guilet et al.[[2015[ ). Because of these complications 
we will not investigate this effect here and refer the reader to the 
studies of [Masada et al.[ ( [2007j[2012[ ) and [Guilet et al.[ ( [20T5] ). 


2.3 MRI termination 

2.3.1 Termination scenarios 

MRI channels cannot grow indefinitely, because their energy would 
constantly increase, whereas the energy of the system contained in 
differential rotation, which provides the reservoir for the MRI, is 
finite. Hence, there must be a physical mechanism terminating MRI 
growth. 

GX94 suggested that MRI channels being exact solutions 
of the MHD equations (in the shearing sheet and incompressible 
flow approximations) may be unstable against parasitic instabili¬ 
ties, which could terminate the MRI growth. They found in their an¬ 
alytic calculations (under the assumptions described in Sec. [2.3.2) 
that in ideal MHD (shear driven) KH modes can develop on top 
of MRI channels. GX94 also suggested that in resistive MHD, par¬ 
asitic instabilities of the (current driven) TM type could develop, 
too. We note that the importance of magnetic reconnection for MRI 
termination was already discussed by BH91. Analytic calculations 
by [Latter et al.[ ( [20()^ in resistive MHD confirmed this hypothesis. 
Alternatively, if this scenario does not hold, the magnetic field of 
the MRI channels could grow to a dynamically relevant strength 
(when the Alfven speed becomes comparable to the sound speed), 
violating the approximation of incompressibility. Consequently, the 
magnetic field pressure would become important and it could push 
matter towards magnetic null surfaces of the MRI channels, or MRI 
channels could become buoyancy unstable (GX94). Pinally, small 
amplitude MRI channels emerging in an already (MRI-driven) tur¬ 
bulent state could be destroyed by non-linear mode-mode interac¬ 
tions or turbulent mixing ( [Latter et al.|2009[ ). 

In this paper, we only consider the first scenario, i.e. MRI ter- 
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mination via parasitic instabilities. In Sec. |2.3.2| we will briefly 
discuss assumptions and flndings of the parasite model of [Good- 1 
man & Xu|{|1994|) and of its extension to resistive-viscous MHD by 
PessahU2010|K 


2.3.2 Termination via parasitic instabilities 

GX94 considered perturbations in a system with already well devel¬ 
oped MRI channels (of the fastest-growing mode with kuRi) given 
by 

V = -qTloir - ro)0 Vc(t; ^mri) + Vp(r, 0, z, t), (33) 

b = bo^z + bc(t; hum) + bp(r, 0, z, t), (34) 

where Vp and bp are the velocity and the magnetic held of the par¬ 
asitic instabilities, respectively. 

Solving the equations governing the evolution of the sec¬ 
ondary (parasitic) instabilities is a very challenging task, because 
MRI channels, which are treated as a background held for the per¬ 
turbations, are non-stationary. Hence, standard techniques like an 
WKB ansatz cannot be used. To make this task more tractable for 
analytic studies, GX94 considered a stage of MRI growth when the 
amplitude of the MRI channels is much larger than the initial weak 
magnetic held, i.e. be » bo^. The growth rate of the secondary in¬ 
stabilities 7 p (which scales oc be) is then much larger than the MRI 
growth rate, i.e. yp » yMRi- Under these conditions the time evo¬ 
lution of the MRI channels, the Coriolis force (which is of the or¬ 
der of yMRi), the background shear flow, and the initial background 
magnetic held bo^ can be neglected. Hence, instead of searching for 
solutions to perturbations according to Eqs. ( [^ and ( [^ , GX94 
considered a more simplifled system where the velocity and the 
magnetic held are given by 

V(0 = Vc(to; burn) + Vp(r, 0, Z, t), (35) 

b(0 = he(to; hum) + bp(r, 0, z, t), (36) 

with to = const, being the time at which the secondary perturbations 
are imposed. Similar assumptions were also made by [Latter et al.| 
( |2009| ) and |Pessah| ( |2010| ). 

( |2010| ) identifled regions in parameter space, where de¬ 
pending on the values of the hydrodynamic and magnetic Reynolds 
numbers either KH or TM is the dominant (i.e. faster developing) 
secondary instability that terminates MRI growth. In particular, for 
the conditions prevailing in core collapse supernovae outside the 
PNS (Re, Rm ^ 1) the exponential growth phase of the MRI should 
be terminated by KH instabilities. TM should be dominant only 
in very resistive media, i.e. if R^ < 1. jPessah] ( |2010| ) also found 
that, in general, (shear driven) KH modes grow fastest along the 
MRI velocity channels, 0 kh = 0v. whereas (current driven) TM 
grow fastest along the magnetic held channels, 0 tm = (pb- One of 
the aims of this work is to test these predictions through numerical 
simulations. 


Pessah 


the code employs various optional high-order reconstruction al¬ 
gorithms including a total-variation-diminishing piecewise-linear 
(TVD-PL) reconstruction of second-order accuracy, a third-, fifth-, 
seventh- and ninth-order monotonicity-preserving (MP3, MP5, 
MP7 and MP9, respectively) scheme ( jSuresh & Huynh] jl 997 j l, 
a fourth-order, weighted, essentially non-oscillatory (WEN04) 
scheme ( [Levy et al.|2002| l, and approximate Riemann solvers based 
on the multi-stage (MUSTA) method ( [Toro & Titarev||2006| ). We 
add terms including viscosity and resistivity to the flux terms in the 
Euler equations and to the electric held in the MHD induction equa¬ 
tion. We treat these terms similarly to the fluxes and electric fields 
of ideal MHD, except for using an arithmetic average instead of an 
approximate Riemann solver to compute the interface fluxes. The 
explicit time integration can be done with Runge-Kutta schemes of 
first, second, third, and fourth order (RKl, RK2, RK3, and RK4), 
respectively. 

Choosing appropriate numerical schemes for our 3D simula¬ 
tions is an important issue, because we want to keep the numerical 
viscosity and resistivity as low as possible without unnecessarily 
increasing the computational cost. Therefore, before applying the 
code to the MRI, we assessed its numerical viscosity and resistiv¬ 
ity in an extended set of auxiliary simulations (Rembiasz et al., in 
preparation). We performed test calculations of linear-wave prop¬ 
agation and the TM instability, comparing the numerical solutions 
with (semi-) analytic solutions. Our tests show that the very low 
numerical dissipation of the MP9 scheme well justifies its larger 
stencil (requiring more ghost zones). In the TM simulations, the 
main contribution to the numerical dissipation comes from the spa¬ 
tial rather than the temporal discretisation errors. We do not ob¬ 
serve any gain in accuracy when using the RK4 instead of the 
RK3 scheme. Because we expect these flndings to hold also for 
MRI simulations, we performed the simulations reported here with 
the MP9 scheme, a MUSTA solver based on the HELD Riemann 
solver, and an RK3 time integrator ( |Harten|1983| [Miyoshi & Ku-j 
|sano|2005] ). 


3.2 Equation of state 

We use the hybrid equation of state (EOS) of jKeil et^ ( |1996| ), in 
which the gas pressure P results from the addition of two contri¬ 
butions, namely a baryonic one Pb and a thermal one Pth- These 
pressure contributions are given by 

= Kp^\ (37) 

/'.h = (r,h - i)eth, (38) 

where K = 4.897 x 10^"^ is the polytropic constant, and Tb = 1.31 
and Tth = 1.5 are the barotropic index and the thermal adiabatic in¬ 
dex, respectively. The quantity ^th is the thermal part of the internal 
energy i.e. ^th = ^ - T’b/(rb - 1). 


3 METHOD 
3.1 Code 

We use the three-dimensional Eulerian MHD code Aenus ( jOber-j 
|gaulinger||2008| l to solve the MHD equations The code 

is based on a flux-conservative, finite-volume formulation of the 
MHD equations and the constrained-transport scheme to main¬ 
tain a divergence-free magnetic held ( [Evans & Hawley] 1988|l. Us - 
ing high-resolution shock-capturing methods (e.g. |LeVeque[1992] |, 


3.3 Computational grid and boundary conditions 

Our study comprises a set of two-dimensional axisymmetric simu¬ 
lations and a set of three-dimensional simulations. For both kinds 
of simulations we employed cylindrical coordinates (r, 0, z) and a 
computational domain centred around the equatorial plane at a ra¬ 
dius ro = 15.5 km. For this value of ro our computational box is 
located in the middle of a nascent PNS of radius rpNs ~ 30 km. 
However, [Guilet et “^ ( [20151 ) have recently shown that the viscos¬ 
ity due to neutrinos can be much higher in collapsing cores than 
previously thought. Close to the neutrinosphere, i.e. whenever the 
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neutrino mean free path is shorter than the length scale of inter¬ 
est, very low Reynolds numbers are expected. The most favuorable 
place for MRI amplification is located close to the PNS star surface, 
where differential rotation is stronger, and the Reynolds numbers 
are larger if the surface is located above the neutrinosphere, which 
however may not always be the case (see Fig. 10 in|Guilet et al.| 
\m5\ . However this recent finding does not invalidate our studies 
as our models can be easily scaled for different initial conditions. 
We may choose different values of the radius, r, rotational velocity, 
Q, or density, p, and scale the key physical quantities in following 
way: 


Mr. = MrJ-l 1^1 1 ^ 


^0 / \^o/ \Po 




- /r \(Q 


1/2 


t = t 


Tmri — Tmri 


Q 


dlVlRI — d^MRI I , 

Vol 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


where Qq = 1824 s"^ and po = 2.47 x 10^^ g cm"^, L/ is the box 
length in the direction i (where i = r, 0, z), and the barred quanti¬ 
ties are the values set up in or computed from our simulations. We 
stress again that the most favorable place for the development of the 
MRI is the region close to the surface of the PNS, where the differ¬ 
ential rotation gradient and the Reynolds numbers are larger than 
deep inside the PNS. Since the choice of Go is ad hoc, it is evident 
that the ratio G/Gq can be made as close as desired to one. This 
means that neither the MRI growth rate nor the typical growth time 
will change under a translation of the box to the surface of the PNS. 
Such a translation will also imply that r/ro ~ 2, while p/po ~ 0.1, 
hence, the typical magnetic fields and Maxwell stresses would be 
0.6 and 0.4 times smaller than in our models. However, as we will 
show, none of these variations would change the foremost qualita¬ 
tive prediction of our models, namely, that the termination of the 
MRI growth in collapsing stellar cores happens by the action of 
parasitic KH modes. 

In the 3D simulations, the typical box size is x x = 
1 km X 4 km x 1 km in r, 0, and z direction, while the typical box 
size is 1 km x 1 km in the (r,z) plane in our 2D simulations. We 
performed simulations involving up to 200 X 800 x 200 zones (see 
Table[T). 

We assume periodic boundary conditions in both 0 and z- 
direction. This choice is natural for the 0 direction, whereas in the 
vertical direction, z, the core is obviously not periodic. However, 
for the simulated region near the equatorial plane the vertical com¬ 
ponent of the gravitational force, Fg^, can be neglected, because it 
is much smaller than the radial one, i.e. Fg^/Fg^ < 0.03. 

Unlike in local simulations of accretion discs, we cannot use 
the shearing sheet boundary condition (see |Hawley & Balbus|1991| ) 
in the radial direction, because it does not allow for global gradi¬ 
ents of thermodynamic variables (which are present in core col¬ 
lapse supernovae) in the simulation domain. Therefore, like |Ober^ 
Igaulinger et al.| ( |200'^ (who followed |Klahr & Bodenheimer|2003| ), 
we use the shearing disc boundary condition, which allows one to 
take these gradients into account. 


We solve the full compressible MHD equations and do not 
perform a transformation to the frame corotating with the fluid, 
which discriminates our simulations from the common shearing- 
box approach. In radial direction, we apply periodic boundary con¬ 
ditions to the deviation of a variable (here, density) from its initial 
background state, e.g. 

6p{r, t) = p{r, t) - p(r, 0), (45) 

i.e. we enforce periodicity of the perturbations. We apply these 
boundary conditions to angular velocity, density, momentum, and 
entropy. Because the initial magnetic field is homogeneous in all 
our simulations, we use periodic boundary conditions for this quan¬ 
tity too. 


3.4 Initial conditions 

Like |Obergaulinger et al.| ( [200^ , we use equilibrium initial models 
based on the final stages of post-bounce cores from|Obergaulinger| 
|ir^ ( [2006^ , in which (several tens of milliseconds after core 
bounce) the shock wave has reached distances of a few hundred 
kilometres and the post-shock region exhibits a series of damped 
oscillations as the PNS relaxes into a nearly hydrostatic configura¬ 
tion. 

The rotational profile (given by Eq.pT] with q - 1.25) that 
we used in our simulations, is similar to the one employed in the 
global MRI simulations of |Obergaulinger et al.| < |200^ . Because 
the resulting centrifugal force is insufficient to balance gravity, the 
gas is kept in (an initial hydrostatic) equilibrium by an additional 
pressure gradient, so that 

- drP + rpQ^ = 0. (46) 


The initial distributions of angular velocity, density, and gravita¬ 
tional potential are depicted in Fig.[2 

Unless otherwise stated, we set the initial magnetic field 
strength to = 4.6 X 10^^ G, the shear and bulk viscosity to 
y = ^ = Ocm^ s“^ and the resistivity to q = 4.45 X 10^ cm^ s“^ 
which implies a Reynolds number Rq - oo and Fm ~ 100, re¬ 
spectively. The Reynolds numbers slightly vary in the box from 
Fm = 89-125, because both the Alfven speed and the rotational 
velocity are functions of radius (see Fig. □ and Eqs.[^ and [^ . 
For these default parameters, the wavelength of the most unsta¬ 
ble MRI mode ranges from Tmri = 0.314-0.385 km, and the cor¬ 
responding velocity channels and the magnetic field channels form 
(for a rotational shear q = 1.25) at an angle, 0v = 44.4-44.6° and 
(pb - 134.7-134.8°, respectively (see Eqs. ( [^ and ([^). These an¬ 
gles differ from those in ideal MHD (0^ = 45° and 0^ = 135°) only 
very little. 

In all simulations presented below the initial magnetic field 
has only a uniform component in z direction, as defined in Eq. 

This field geometry is a popular choice in MRI simulations, (see, 
e.g. BH91;|Hawley & Balbus|1991|[Sano & Inutsuka|2001[|Ober-| 


Igaulinger et al.||20'09 l, since the vertical component is the most 

important one for the development of the instability (cf. |Balbus| 
|& Hawley|1998| l. Another common choice is a so-called zero net 
flux configuration (see, e.g. Fromang et al. 12007) |Fromang & Pa-| 


|paloizou|2007[ |Obergaulinger et al.|2009j l, in which the magnetic 

field has a sinusoidal radial dependence, i.e. b oc z sin(^^r), where 
kr is chosen in such a way that an integer number of wavelengths 
fits the computational domain. Thus, K - InnlLr with n being a 
natural number. 

The value of the initial magnetic field amplitude, requires 
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Figure 1. Hydrostatic structure of the initial models. The diagram shows 
the gravitational potential O 19 = 0 /( 10 ^^ erg cm“^) (solid line, right ordi¬ 
nate), the density pi 3 = p/( 10 ^^ g cm“^) (dashed line, left ordinate), and 
the angular velocity ^3 = n/(10^ s“^) (dash-dotted line, left ordinate). The 
entropy profile of this specific model is assumed to be fiat. 


We pay attention to choose the values of the initial magnetic 
field strength such that an integer number of channels of the fastest- 
growing MRI mode fits into the computational domain. In this way 
this mode fulfils the periodic boundary conditions imposed in z di¬ 
rection. Otherwise, the fastest-growing mode could be artificially 
suppressed by an unfavourable box size. |Rembiasz| pOTS] ) showed 
in test simulations that if, e.g. a box of size 2.5 /Imri in z direction 
is used, usually three MRI channels form. Either all three of them 
have a wavelength smaller than Tmri or he finds a combination of 
larger and smaller channels. In any case, the MRI developed at a 
rate which was lower than theoretically expected. 

To trigger the MRI we impose an initial velocity perturbation 
on the background velocity profile (defined by Eq.|^ of the form 

Vi = Qr[{{b^9t^(r, 0,z) + 6 sin(^^z)}f}-i- 

{6%{r, 0 , z)}0 + 0 , z)z], (48) 


where 9t^(r, 0 , z), 0 , z), and 0 , z) are random numbers in 

the range [-1,1], S and Sr are the perturbation amplitudes, is 
the radial perturbation wavenumber, and e is the amplitude of the 
sinusoidal perturbation. If not ot herwise written, - kum, Sr = 
10"^, S = 10"^, and e = 2x 10"^. Obergaulinger et al. ( 2009) used 
a similar prescription for the initial perturbation, except that the 
sinusoidal part was not present in their case, i.e. 6 = 0. We find 
that the sinusoidal term is more robust in exciting MRI modes from 
small perturbations, which are sometimes suppressed by numerical 
effects if only random perturbations are imposed (see |Rembiasz| 
|20T3l ). 


some further comments. According to state-of-the-art stellar evo¬ 
lution calculations the pre-collapse magnetic field for the most 


strongly magnetised progenitors is less than about 10^ G (Heger 


|et al |20()5] l. During the collapse phase the magnetic field can be 
amplified by compression by two orders of magnitude to ^ 10^^ G 
( [Meier et al.|I976| l. Erom Eq. ( [^ , which is valid in ideal MHD, 
we estimate that for the PNS the wavelength of the fastest-growing 
MRI mode is 


A 


MRI 


^ 70 cm 


'l / P 

10i‘GjU-5xl0‘3gcm-3 



n r‘ 

1900 s-i ) ■ 


(47) 

Eor the typical Reynolds numbers used in our simulations (Re = 
Rm ~ 100) the wavelength of the fastest-growing mode is ^ 1.5% 
longer (and the growth rate ^ 1 . 6 % lower) than the one in ideal 
MHD. [^Assuming that 10 zones per MRI channel are needed to 
resolve it properly, a simulation with boz = 10^^ G would require a 
resolution of the order of 10 ^ zones per dimension, which is already 
unaffordable in 2D. One could reduce the cost by using a smaller 
computational domain, but the high rotational and sound speeds 
(v 0 ^ Cs ~ 3 X 10^ cm s"^) would still limit the timestep to At ^ 
3x10“^ ms, i.e. almost 10 ^ iterations would be required to simulate 
the MRI until its termination at ^ 15 ms. 

Therefore, following [Obergaulinger et aL] ( |2009j ), we use ini¬ 
tial magnetic fields which are two orders of magnitude higher 
(boz ~ 10^^ G) to increase the wavelength of the fastest-growing 
MRI mode to ^ 100 m. This reduces the minimum resolution to 
^ 100 X 400 X 100 zones in a 1 km x 4 km x 1 km box, and the 
number of iterations to less than 10 ^ . 


^ We obtained these values by plotting 7 mri(^mri) (for given Q,q,CAz,^, and 
y) using the expression from PC08 and determining the maximum 7 mri and 
its location, kMRi, graphically. 


4 RESULTS 


4.1 2D simulations 

4.7.7 Termination in 2D 


Imposing axisymmetry severely limits the number of modes that 
can grow in 2D simulations. While the fastest-growing MRI mode, 
which is an axisymmetric one, can freely develop in such sim¬ 
ulations, the dominant parasitic instabilities, which for Rm > 1 
are non-axisymmetric KH modes, are suppressed ( [Pessah[[2010[ l. 
Hence, among axisymmetric secondary instabilities the fastest- 
growing modes are of TM rather than KH type. This even holds 
for simulations with a very low or even a vanishing physical resis¬ 
tivity, as discussed by ( [Obergaulinger et al.|2009[ ). These authors 
performed extensive studies of the MRI by means of local 2D ideal 
MHD simulations. Their simulations confirmed the instability cri¬ 
teria and the growth rates of the MRI for the flow regimes relevant 
to core collapse supernovae. They also found that the growth of the 
MRI is terminated by a tearing mode (TM) instability developing 
because of the unavoidable presence of a numerical resistivity in 
(even ideal) finite-volume MHD codes. 

Pigure|^ summarizes the evolution of an axisymmetric model 
simulated with a resolution of Nr = = 100 zones. We performed 

appropriate convergence studies to ensure that the MRI is properly 
resolved at this grid resolution (see [Rembiasz|[2013[ for details). 
The top left panel displays the time evolution of the absolute value 
of the volume-averaged Maxwell stress component 


-NXrcf) — 


1/ brb^ dv\ 

V 


(49) 


where V is the volume of the computational domain. The other 
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Figure 2. Axisymmetric (2D) MRI simulation in which the instability is terminated by TMs. Top left: Time evolution of the absolute value of the volume 
averaged Maxwell stress component Atr^ Eq. j^ . The three green vertical lines mark the times corresponding to the snapshots shown in the other three 
panels, which display the structure of the radial magnetic field at f = 4ms (top right), f = 12.1 ms (bottom left), and f = 13 ms (bottom right), respectively. 


three panels display the colour encoded value of the radial com¬ 
ponent of the magnetic field in the r - z plane at three different 
times. 

From the initial velocity perturbations (see Eq.[^, three MRI 
channels have formed at ^ = 4 ms {top right panel) that grow expo¬ 
nentially with time at a constant rate {top left panel). A linear fit of 
\ogMr(p{t) in the time interval t G [6,8 ] ms gives yMRi = 1127 s"^ 
This value is consistent with the local linear analysis, Eq. ( [25] ), 
which predicts values of 7 mri varying from 1087 s“^ to 1175s"^ 
within the box boundaries because of its dependence on i2(r). 

Ait - 12 ms, the Maxwell stress reaches a value of Mrcp - 


2.78 X 10^^ and the MRI growth is terminated by parasitic in¬ 
stabilities. We note that in 2D simulations, the stress at termi¬ 
nation is highly sensitive to the initial random perturbation im¬ 
posed in the simulation. Performing several realizations of the same 
simulation, with different seeds for the random number genera¬ 
tor, we obtained a Maxwell stress at the termination varying from 
Mrcp = 2.19 X 10^0 to Mr^ = 4.69 X 10^0 G^ (see 

for details). As we discuss in the next section, we do not observe 
this large scatter in our 3D simulations. 

The colour map of the radial component of the magnetic 
field exhibits several X points at ? = 12.1 ms {bottom left panel). 


Rembiasz 2013 
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Figure 3. Relation between the (not necessarily orthogonal) coordinates 
ih,w) defined by Eqs. and (HI and the cylindrical coordinates (r,z). 
The angles and 0^ are given by Eqs. and respectively. Compare 
withFig.2ofPC08. 


where field lines reconnect. These X points are located at (r, z) = 
(15.2,0.4), (16,0.1), and (15.2,-0.3) km, respectively. There are 
also three O points recognisable in the centres of magnetic islands 
located at (r,z) = (15.5,0.4), (15.7,0.1), and (15.8,-0.3) km, re¬ 
spectively. These patterns indicate MRI termination by TM rather 
than by KH instabilities. Shortly afterwards, at ^ = 13 ms {bottom 
right panel), the channel modes have been destroyed and MHD tur¬ 
bulence sets in. 

These results, which are qualitatively similar to those found 
by |Obergaulinger et al.| ( |200^ , confirm that in axisymmetric mod¬ 
els the MRI is artificially terminated by TMs, which would grow 
slower than KH instabilities in 3D models. 


4.2 3D simulations 

4.2.1 Termination in 3D 

According to predictions of |Pessah| ( |2010| ), in general, the fastest- 
growing KH instabilities should develop along the velocity chan¬ 
nels, i.e. 0 kh = 0v. whereas the fastest TM should grow along 
the magnetic field channels, i.e. 0 tm = 0 ^ 7 - As these channels de¬ 
fine two important directions in the horizontal (r, 0 ) plane, we will 
sometimes use in the following discussion another coordinate sys¬ 
tem (/t, w, z), the axes h and w being aligned with the velocity and 
magnetic field channels, respectively. The {h,w) coordinates and 
the (r, 0 ) cylindrical coordinates are related through the following 
coordinate transformation (Fig.[^ 

^ = rcos0v + 0sin0v, (50) 

= rcos0^ -h 0sin0^, (51) 

where the angles 0v and 0 ^ are given by Eqs. ( [^ and ( [^ , respec¬ 
tively!^ For Re = the axes h and are orthogonal, i.e. during 
the phase of exponential growth of the MRI bh-0 and = 0. 

^ This naming convention differs from that of |Latter et al.H2009t and |Pes-| 
|sah| ( |201()) who used h to denote vectors (or their components) in the (r, 0) 
plane. 



Figure 4. Sketch of the computational domain of size 1x4x1 km, which 
was used to simulate models #5 to #8 and #10 (see Table[T} together with 
two cuts corresponding to the {h, z) plane (blue) and the (w, z) plane (red) 
used to display some results of model #7 in Fig.[^ Note that the computa¬ 
tional volume is actually bent in 0 direction and the two cuts are no planes 
but curved surfaces in cylindrical coordinates. 


Figure!^ depicts the horizontal plane in both coordinate sys¬ 
tems. The channel angles are set to 0v = 45° and 0 ^ = 135°, which 
corresponds to the ideal MHD limit. In general, the vectors h and 
w are not orthogonal to each other, and their exact orientation de¬ 
pends on the Reynolds numbers. For example, in the limit ^ oo 
and Rm ^ 0, one has h\\r and h'H^, i.e. 0v = 0° and 0 ^ = 90° 
(PC08). 

For the values of the Reynolds numbers used in our simu¬ 
lations, i.e. Re = oo,Rm ~ 100 , and Re = Rm ~ 100 , we ex¬ 
pect the channels to be oriented along 0v = 44.4-44.6°, 0 ^ = 
134.7-134.8°, 0v = 44.3-44.5°, and 0 ^ = 134.5°, respectively. 
These angles differ only little from those of the ideal MHD case. 
Hence, in either case, the fastest-growing KH mode should develop 
at an angle 0 p = 0 kh ~ 44.5°, where 0 p denotes the angle at which 
the dominant parasitic instability develops. Depending on the initial 
conditions, one either has 0p = 0 kh or 0 p = 0 tm- In 2D axisymmet¬ 
ric simulations, the only allowed angle is 0 p = 0 °. 

We studied the termination process in a number of 3D simu¬ 
lations, varying the size and aspect ratio of the computational do¬ 
main, and the grid resolution (see Tab.[^for the list of models). The 
default simulation box has a size of L^xL^xL^ = 1 kmx4 kmx 1 km, 
as shown in Fig.j^ which also depicts two cuts corresponding to 
the (/t, z) and (iv, z) planes. We note that the computational volume 
is actually bent in 0 direction and the two cuts are not planes but 
curved surfaces in cylindrical coordinates. However, we will call 
them planes for simplicity. 

We begin the discussion of our results with model #7 (see 
Tab.[T}, which we simulated in the default box resolved with 
100 X 400 X 100 zones in r, 0, and z direction, respectively. Except 
for the third space dimension, the initial conditions and the param¬ 
eters of this 3D model are identical to those of the axisymmetric 
model discussed in Sec. |4.1| As illustrated by the time evolution of 
the Maxwell stress (Eig.|^, the MRI grows exponentially with time 
at the same rate as in the 2D model ( 7 mri = 1127 s"^) until satura¬ 
tion, which occurs a bit earlier than in the 2D model att - 11.2 ms. 

Eigure|^ shows the 3D structure of the radial component of 
the magnetic field at six distinct times (marked by vertical lines in 
Fig.|5). Magnetic field perturbations grow in the form of three ax¬ 
isymmetric channels {upper left panel), which are perturbed in turn 
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Table 1. Overview of our 3D MRI simulations, which were all performed with the initial parameters chosen in such a way that /Imri ~ 0.333 km*. The 
columns give the model identifier, the magnetic field strength, the hydrodynamic and magnetic Reynolds numbers (Re and Rm, respectively), * the size of the 
computational domain, the resolution, the number of grid cells per MRI wavelength, the measured MRI growth rate, the volume-averaged Maxwell stress at 
termination, the azimuth angle at which the parasitic instability develops, and its horizontal wavelength (i.e. in the (r, 0) plane) measured at f = 11 ms and 
t = 12.5 ms for KH instabilities and TMs, respectively. The azimuth angle and the horizontal wavelength are determined with about a 9 ° and 10% accuracy, 
respectively (see also Fig. EH. The final column gives the type of the instability. For models 19 and 20, we were unable to identify the type of the parasitic 
instability terminating the growth of the MRI. 


# 

boz 

[10‘3 G] 

Rc 


box size 

(r X 0 X z) [km] 

resolution 
(r X 0 X z) 

zones per 

channel 

7mri[s M 

AXrtp 

[10^® G^] 

0p [°] 

A 

^MRI 

term. 

instab. 

1 

4.6 

oo 

oo 

1x4x1 

100 X 400 X 100 

33 

1137 

0.96 

45 

0.94 

KH 

2 

4.6 

oo 

100 

1 X 1 X 0.333 

24 X 24 X 8 

8 

1104 

0.79 

45 

1.5 

KH 

3 

4.6 

oo 

100 

1 X 1 X 0.333 

30 X 30 X 10 

10 

1122 

1.3 

44 

1.2 

KH 

4 

4.6 

oo 

100 

1 X 1 X 0.333 

48 X 48 X 16 

16 

1130 

1.1 

47 

1.1 

KH 

5 

4.6 

oo 

100 

1x4x1 

60 X 240 X 60 

20 

1126 

1.1 

44 

0.92 

KH 

6 

4.6 

oo 

100 

1x4x1 

76 X 304 X 76 

25 

1127 

1.0 

47 

0.99 

KH 

7 

4.6 

oo 

100 

1x4x1 

100 X 400 X 100 

33 

1127 

0.93 

44 

0.85 

KH 

8 

4.6 

oo 

100 

1 x1x1 

100 X 100 X 100 

33 

1127 

0.93 

47 

0.95 

KH 

9 

4.6 

oo 

100 

1 X 1 X 0.333 

100 X 100 X 34 

34 

1127 

0.93 

47 

0.77 

KH 

10 

4.6 

oo 

100 

1x4x1 

200 X 800 X 200 

67 

1127 

0.73 

44 

0.70 

KH 

11 

4.6 

oo 

100 

1 X 1 X 0.333 

400 X 400 X 134 

134 

1128 

0.73 

45 

0.69 

KH 

12 

4.6 

100 

100 

1 X 1 X 0.333 

100 X 100 X 34 

34 

1120 

0.90 

44 

0.85 

KH 

13 

4.6 

100 

100 

1 X 0.8x0.333 

100 X 80 X 34 

34 

1120 

0.92 

44 

0.80 

KH 

14 

4.6 

100 

100 

1 X 0.6x0.333 

100 X 60 X 34 

34 

1120 

1.0 

43 

0.96 

KH 

15 

4.6 

100 

100 

1 X 0.4x0.333 

100 X 40 X 34 

34 

1120 

0.97 

43 

0.80 

KH 

16 

4.6 

100 

100 

1 X 0.3 X 0.333 

100 X 30 X 34 

34 

1120 

1.1 

52 

0.70 

KH 

17 

4.6 

100 

100 

1 X 0.2x0.333 

100 X 20 X 34 

34 

1120 

4.5 

- 

- 

TM 

18 

4.6 

100 

100 

1 X 0.1 X 0.333 

100 X 10 X 34 

34 

1120 

4.9 

- 

- 

TM 

19 

0.325 

oo 

0.1 

1 X 1 X 0.333 

100 X 100 X 34 

34 

80 

0.0020 

- 

- 

7 

20 

0.163 

oo 

0.05 

1 X 1 X 0.333 

100 X 100 X 34 

34 

34 

0.0011 

- 

- 
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* Note that Tmru Rq, and Rm are not uniform throughout the computational domain, but vary by ^ 20% (see Sec. |3.4| for details). 


by secondary instabilities that become recognisable only shortly 
before the MRI is terminated at t = 11.2 ms. The sequence of snap¬ 
shots around this time {upper middle to bottom middle) demon¬ 
strates that non-axisymmetric parasitic modes are responsible for 
MRI termination. After exponential growth of the magnetic field 
has ended, the simulation volume is dominated by a small-scale, 
turbulent magnetic field {bottom right). Comparing the structure of 
the radial magnetic field during MRI termination (Fig.|^, we find 
profound differences between the 2D model {left-hand panel) and 
the 3D model {right-hand panel). The 2D model exhibits X and O 
points that are characteristic of the TM instability, whereas in the 
3D simulation the channel modes bend strongly at the locations of 
vortex rolls which are a typical feature of the KH instability. 

We also studied some geometrical aspects of the parasitic in¬ 
stabilities found in the 3D simulation. For this purpose, we use the 
(/t, w, z) coordinates.Velocity-shear driven KH instabilities should 
grow fastest in the h direction and dominantly current driven TM 
along the axis. To simplify the expressions we performed the 
coordinate transformation (r, 0 , z) {h, w, z) with 0 v = 45° and 
0 ^ = 135° instead of the theoretically expected angles 0v = 44.5° 
and 0^ = 134.8°, which does not affect our qualitative analysis, 
however. 

Figurej^ shows the distribution of the magnetic field com¬ 
ponents bh {left-hand panels) and b^ {right-hand panels) in two- 
dimensional cuts through the computational domain of the 3D 
model#? shortly before MRI termination. The component bh is 
considerably smaller than the component b^. This is consistent 


with the theoretical expectation for an MRI channel whose mag¬ 
netic field should grow in the iv direction (0^ = 135°) and vanish in 
the perpendicular h direction (0 = 45°). Close to termination, the 
MRI channels are strongly perturbed and vortex rolls begin to form 
in the h direction {upper panels). This indicates that KH instabili¬ 
ties developing along the velocity channels (i.e. in the h direction) 
are responsible for the disruption of the MRI channels. In the w di¬ 
rection (i.e. along the magnetic field channels separated by current 
sheets in which TM could develop), the channels remain almost un¬ 
perturbed. Only box-size structures appear which are most likely a 
consequence of the radial boundary conditions imposed in our sim¬ 
ulations (see Sec. |3.3| ). Although we cannot discard the presence of 
sub-dominant TM in this projection, it is clear that TM do not play a 
dominant role in the termination process, which can be understood 
completely in terms of parasitic KH instabilities. 

Before we proceed further we take another look at the spatial 
structure of br during MRI termination at ^ = 11.1ms, which is 
illustrated in Fig.|^ In the upper panel depicting br in a (0, z) cut 
at r = 16 km (this cut corresponds to the outer radial boundary of 
the computational domain shown in the upper right panel of Fig.[^, 
one can recognize vortex rolls developing along the MRI channels. 
From the middle panel, which shows the distribution of br in the 
(r - 0) plane at z = 0.335 km, one could determine the horizontal 
wavelength Tp and the angle at which parasitic instabilities develop 
(0p 45°). However, we can obtain these quantities much more 

accurately using Fourier transforms (see next subsection). The bot¬ 
tom panel provides another view of the structure of br. By showing 
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time [ms] 

Figure 5. Evolution of the volume-averaged Maxwell stress component 
Mr(p for the 3D models #5 (orange), #6 (red), #7 (black), #10 (violet), and 
#11 (blue; run only until t = 12 ms). Up to r ^ 10.5 ms, the models give 
almost identical results, differences being only visible around termination 
(see inset). Vertical lines indicate the times of the snapshots displayed for 
model#? in Fig.[^ For a better readability of the figure, the vertical lines 
around MRI termination are only displayed in the inset. 


a part of both the (h, z) cut and the (w, z) cut, it illustrates in which 
direction parasitic instabilities develop. We note that the (h, z) and 
(w, z) cuts shown in Figs. andare different. 

4.2.2 Fourier analysis 

In order to confront the theoretical expectations with our numer¬ 
ical results, we have calculated spatial discrete 3D Fourier trans¬ 
forms of the magnetic field components ba with cr G r, 0 at a given 
time using a fast Fourier Transform (FFT) algorithm. We denote the 
complex FFT coefficients as The power spectral density is pro¬ 
portional to \aa\^, which is a measure of the average magnetic field 
energy density of the component ba in Fourier space. We chose 
the components br and b^p for the analysis, because they contain 
the information about both the MRI channel flows and the para¬ 
sitic instabilities. We are particularly interested in determining the 
wavevectors, k = {kr^k^p^kz), of the dominant modes, which have 
the largest amplitudes in the Fourier space. We expect that MRI 
channel flows appear as structures with a wavevector 

kmri = (0, 0, ka^i), (52) 

whose modulus is maximum for the fastest-growing mode, ^mri- 
Parasitic instabilities develop in the whole Fourier space, 
whereas the MRI does not contribute to modes with finite kr and 
ktp, but kz - 0. Hence, parasitic instabilities are expected to produce 
a characteristic signature with wavevectors 

kp = (^^,/:^,0), (53) 

which should be distinguishable from the MRI signature. The max¬ 
imum Fourier amplitude should be attained for the wavevector of 
the fastest-growing parasitic mode. 


The analysis of model#? reveals that the Fourier amplitude 
along the line kr - k^p - 0 peaks at kz - 18.8 km"^ which 
corresponds to the wavelength of the fastest-growing MRI mode, 
^mri ~ 0.333 km of this model. The latter result holds for all our 
simulations during the phase of exponential growth of the MRI. 
The Fourier amplitudes in the plane kz - 0, displayed for three dif¬ 
ferent times for model#? in Fig.[^ show a power excess which 
peaks at {K^k^p) ^ (0.?, 0.?) km The location of this maximum 
barely changes with time, while the amplitudes of the Fourier coef¬ 
ficients increase relative to those at (0,0, kum)- This result is con¬ 
sistent with a super-exponential growth of parasitic instabilities. In 
addition, as we will demonstrate below, the behaviour of the Fourier 
amplitude in the kz-0 plane is consistent with the development of 
parasitic KH instabilities feeding off the MRI channels. 

The average magnetic energy density of the field component 
ba can be computed from the Fourier amplitudes as 

^ Nr 12 Npl2 NJ2 

^mag,Qf “ 2 ^ ^ ^ \^a(hh^m->bn)\ ■> (54) 

^ 1=-Nrl2 m=-Npl2 n=-Nzl2 

where 


(ki,km,kn) 


ilnl 2nm 2nn\ 
\ Lr L(P Lz I 


(55) 


Similarly, we can estimate the average magnetic energy density of 
the MRI channels or parasitic instabilities restricting the summa¬ 
tion to locations in Fourier space relevant for each kind of instabil¬ 
ity. Therefore, we define 

= ba(0, 0 ,/:mri)P, (56) 

^ Nr /2 Np /2 

ep,a = 2 X X (57) 

^ l=-Nr/2 m=-Npl2 

as proxies for the average magnetic field energy density associated 
with these instabilities. In case of the MRI, the Fourier amplitudes 
are distributed along the line kr - k^p = 0, i.e. should be a 

good estimator. However, in the case of the parasitic modes, e^^a 
contains not all contributions to the energy density and thus pro¬ 
vides only a lower bound of the energy density associated with the 
parasitic instabilities. 

Fig.[^ shows the time evolution of the magnetic energy den¬ 
sity for both the MRI and the parasitic instabilities. The behaviour 
of ^MRi,a is very similar to the time evolution of the Maxwell stress 
(see Fig.|^. It is characterized by an exponential growth with the 
same growth rate as for the Maxwell stress Mrcp, and a termination 
point associated with the disruption of the MRI channel flows. The 
average magnetic energy density e^^a of the parasites starts to grow 
super-exponentially at ^ 9 ms from a value of about 8 orders of 

magnitude smaller than that of the MRI. At termination, the mag¬ 
netic energy density of the parasites amounts to, at least, 6% of the 
magnetic energy density of the MRI. 

To determine the wavevector of the fastest-growing parasitic 
instability, we search for the maximum Fourier amplitude with 
kz = 0. For this purpose, it is sufficient to consider only posi¬ 
tive components of the wavevector k. Numerically, we obtain the 
wavevector by computing the energy-weighted barycenter in the 
relevant part of Fourier space, 

T _ ^l=\^m=\\^r(bhkm,0)\ ka 

“ 2/=li:.=l \ar(kl,km,0r ’ ^ ^ 

where we limited all summations to those Fourier modes displayed 
in Fig.[^to avoid high frequency contaminations. Substituting 
for a^p does not change the results significantly. 













































z [Ian] 


12 T. Rembiasz et al. 







Figure 6. Distribution of the radial component of the magnetic field br across 
times. The times of the snapshot are marked by green vertical lines in Fig.[^ 
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Figure 7. Comparison of the distribution of the radial magnetic field br around MRI termination in a 2D (left) and 3D (right) simulation (model#?; cut at 
(f) = -2 km). Note that the two snapshots are taken at different times. 
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Figure 8. The panels show the distribution of the magnetic field components bh (left) and (right) in two-dimensional cuts through the computational domain 
for the 3D model#? at time t = 11.1 ms. The cut directions are perpendicular (top panels; {h,z) cut) and parallel (bottom panels; (w, z) cut) to the direction of 
the magnetic MRI channels. The locations of the cuts are marked in Fig.j^in blue and red, respectively. 


The values of obtained from Eq. ( [58] ) (white squares in 
Fig.[^ properly trace the location of the maximum Fourier ampli¬ 
tude. With these values one can compute both the wavelength and 
the angle of the parasitic instability according to 

ip = -7^=^ (59) 

= arctan j. (60) 

Figure[T^ shows the time evolution of the wavelength, Tp, and an¬ 
gle, (pp, of the parasitic instabilities during the late stage of MRI 
evolution. Before MRI termination, the evolution of the angle is 
compatible with its theoretically expected behaviour for KH insta¬ 
bilities, i.e. 0p ^ 45° (horizontal green line), within the accuracy of 
the angle determination]^ The wavelength differs from its theoreti¬ 
cally expected value Tp ^ 0.56 km (horizontal blue line) by a factor 
of ~ 2. Whether the source of this disagreement is of a numerical 
origin or results from a limitation of the theoretical approach of 
|Pessah| ( |20I0| ) is beyond the scope of this work. 

^ The accuracy in the determination of the angle and the wavelength, ^ 9° 
and ^ 10% respectively, depends on the accuracy in the determination of 
kp, which is set by the size of the box. 


We have performed a similar analysis for all our other sim¬ 
ulations. Table[^ gives the values of Tp and (pp at a representative 
time before termination (i.e. at 11 ms for models #1 to #16, and at 
12.5 ms for models #17 and #18). In the following three subsec¬ 
tions, we discuss the influence of different numerical and physical 
parameters on the values of these two quantities. 

4.2.3 Box size 

Next, we address whether the size and aspect ratio of the computa¬ 
tional domain influences MRI termination (for a somewhat sim¬ 
ilar study of the post-termination phase, see |Bodo et ^|2008| ). 
Following lObergaulinger et aTT] ( |2009| ), we simulated models #5, 
#6, #7, and #10 in a box of (the default) size (L^ x x = 
1 km X 4 km x 1 km). We performed additional simulations reduc¬ 
ing the size of the box in 0 direction (1 km X 1 km x 1 km, model 
#8), and both in (p and z direction (1 km x 1 km x 0.333 km, models 
#2, #3, #4, #9, and #11). Finally, we computed several models (#12 
to #18) where we varied the azimuthal size L^p of the domain (see 
Table[T}. 

The main motivation for using a smaller box in some of our 
simulations was computational cost reduction. In accordance with 
theoretical predictions of |Pessah| ( |2010| ) for flows with » F 

we found that in model #7 parasitic instabilities develop at an an- 
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Figure 9. Structure of the radial magnetic field component br of the 3D model #7 at f = 11.1 ms. Top: (z, 0) cut at r = 16 km. Middle: (r, 0) cut at z = -0.09 km. 
Bottom: 3D view of the surface of the computational domain where some part of it has been removed to show pieces of the ih,z) and (w,z) cuts. Note that 
these cuts differ from those displayed in Fig.|^ 
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logjo(Normalised Amplitude) at t = 10.9 [ms] 
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Figure 10. Fourier amplitudes of the radial magnetic field br for model #7 
at 10.9 ms (top), 11.1 ms (middle), and 11.2 ms (bottom), respectively. The 
colour coded quantity is the ratio \a(kr, kcp, 0)|/|fl(0,0, ^mri)I of the Fourier 
amplitude representing parasitic instabilities and the one representing MRI 
channels. The white square is the location of the energy-weighted barycen¬ 
ter computed with Eq. 
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Figure 11. Evolution of the average magnetic energy density associated 
with the MRI channels (eMRi,a) and the parasitic instabilities (e^^a) for dif¬ 
ferent components ba of the magnetic field for model #7. 


gle 0p ^ 45° (see Fig.[^ and middle panel of Fig.|^. This result 
suggests that it is unnecessary to use a box elongated in azimuthal 
direction (L^ > L^). Instead, one could rather study the MRI in a 
box with a horizontal aspect ratio to minimize the vol¬ 

ume of the computational domain without affecting the develop¬ 
ment of parasitic instabilities. To test this hypothesis, we simulated 
model #8 in a box of size x x = 1 km x 1 km x 1 km 
using the same spatial resolution as for reference model#?. In both 
simulations, the Maxwell stress tensor at MRI termination attained 
the same value, which confirms our expectations. 

We can further reduce the computational domain in the ver¬ 
tical direction based on the following observation. In simula¬ 
tions performed with the default box, we chose the magnetic field 
strength in such a way that three of the fastest-growing MRI modes 
fit in the computational domain. Thus, it should be possible to 
reduce the vertical extent of the box by a factor of 3, i.e. = 
^mri ~ 0.333 km, without hindering the growth of the dominant 
MRI mode. However, in such a smaller domain parasitic instabil¬ 
ities are restricted to modes having a vertical wavelength equal 
to the width of an MRI channel (called Type-I modes by GX94), 
and modes with a longer vertical wavelength (Type-IF) are sup¬ 
pressed. Nonetheless, according to the predictions of GX94 and 
|Pessah| ( [2010| ), in the ideal MHD limit or for R^^Rm > 1, the dom¬ 
inant parasitic modes should be (KH modes) of Type-I. Hence, we 
expect the MRI termination process to be unaltered by a box of 
smaller vertical size. 

To test this theoretical prediction, we simulated model #9 in a 
box of size x x = 1 km x 1 km x 0.333 km using the same 
grid resolution as in models #7 and #8. We found that the Maxwell 
stress tensor reached the same value at MRI termination in all three 
models. This result not only confirms (given our initial conditions) 
the termination of the MRI by “Type-I” parasitic modes, which is 
an important result, but it also justifies the use of a computational 
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Figure 12. Characteristics of the parasitic instabilities for model #7 during 
the late stage of MRI evolution. Green crosses and blue diamonds depict the 
angle, 0p, and the wavelength, 4, of the parasites, respectively. Horizontal 
solid lines show the corresponding theoretical values for KH modes that 
one expects to develop in this simulation. 


box with a 12 times smaller volume for studying the MRI termina¬ 
tion process. We made use of this fact in our resolution studies (see 
next subsection). 

In Section [TT] we showed that the MRI is terminated by TMs 
in axisymmetric 2D simulations. One can also impose axisymme- 
try in 3D simulations by choosing a box of vanishing azimuthal 
length, i.e. ^ 0. Therefore, we expect that for a certain small, 
but non-zero L^, KH instabilities should be suppressed even in 3D 
simulations, and the MRI should be terminated by TM instead. To 
determine this critical azimuthal box length, we performed a se¬ 
ries of simulations varying from 1.0 to 0.1 km (models #12 to 
#18). Moreover, to investigate the influence of viscosity, we ran 
these simulations with a non-zero shear viscosity corresponding to 
a hydrodynamic Reynolds number Rq - 100. The shear viscosity 
should change the wavelength of the fastest growing MRI mode, 
but only by a negligible amount from the numerical point of view 
(i.e. less than 1%). 

As expected, the MRI is terminated by KH instabilities in 
model #12, which only differs from model #9 by the value of R^. In 
both simulations, the parasitic instabilities have the same horizontal 
wavelength, dp, and they develop at the same angle 0p within the 
measurement error. In models #13-#16, the MRI is also terminated 
by KH instabilities, but the value of the Maxwell stress tensor at ter¬ 
mination is somewhat larger. This result can be easily understood: 
parasitic modes that would grow fastest in a box with = 1 km are 
affected (or even suppressed) in narrower boxes, i.e. the MRI can 
operate a bit longer before it is Anally terminated. The smaller the 
box, the stronger the suppression of the KH modes and the larger 
the deviation of such quantities as the Maxwell stress at termina¬ 
tion from model #12. Finally, in models #17 (L^ = 0.2 km) and 
#18 (L^ = 0.1km), KH modes are suppressed and the MRI is 
terminated by TMs, showing similar features as in the 2D simu¬ 
lations described in section 14.11 We note that the Maxwell stress 
reaches considerably higher values at termination in these models, 
because TM are growing more slowly. A similar behaviour was ob¬ 


served by [Lesaffre et al.| ( |2009| ) in their mean-field shearing-box 
simulations of accretion discs. In simulations done in a box of size 
Lr - L(p - - Tmri, KH modes were suppressed and the MRI 

channels were disrupted by TM, whereas in simulations done in 
larger boxes in the horizontal directions (i.e. Lr - - /Imri and 

L(p - 4/Imri, or Lz - Tmri and Ly - - 4/1mri), the MRI was 

terminated by KH instabilities. 

To summarize, MRI termination in models #12—#18 is well 
explained within the parasite model, and the minimum azimuthal 
box length for which KH modes are not severely affected by bound¬ 
ary conditions is ^ 0.3 km ^ I/Imri- 


4.2.4 Influence of grid resolution 

To properly resolve MRI termination, one needs to use a resolution 
that is high enough not only to resolve the MRI channels, but also 
parasitic instabilities. The latter criterion is more stringent, as par¬ 
asitic instabilities develop finer structures in vertical direction than 
MRI channels do (GX94, |Pessah|2010> . 

We explore this issue in a series of simulations performed in 
the small {Ly x x = 1 km x 1 km x 0.333 km; models #2, #3, 
#4, #9, and #11) and in the default {Ly xL^xL^- I km x 4 km x 
1 km; models #5, #6, #7, and #10) computational box. As we have 
demonstrated in the previous subsection, parasitic instabilities can 
develop equally well in both computational domains, i.e. the boxes 
are equally well suited for resolution studies. 

In models #2—#4, which were simulated with 8-16 zones per 
MRI channel in the vertical direction, we observe differences al¬ 
ready during the phase of exponential growth. First, channel modes 
emerge ^ 2 ms later than in the other better resolved simulations. 
Secondly, the channel modes somewhat differ from the analytic so¬ 
lution, i.e. some imperfections develop on channels which are of 
numerical and not physical origin. Thirdly, the MRI growth rate 
measured with the help of Mycp slightly fluctuates during the phase 
of exponential growth. These fluctuations of the order of a few per¬ 
cent can be explained by the fact that the above mentioned imper¬ 
fections also contribute to My(p. 

In models #5-#ll, which were simulated with at least 20 
zones per MRI channel, the channel modes appear roughly at the 
same time and grow at the same constant rate during the phase 
of exponential growth. Based on these observations, we conclude 
that 20 zones per MRI channel are sufficient to resolve the phase 
of exponential growth with the MP9 scheme (but more zones will 
be required for lower order reconstruction schemes). This result is 
consistent with previous estimations of the numerical viscosity and 
resistivity of the code done by |Rembiasz| ( |2013T ). With the help of 
tests involving Alfen wave propagation and TM instabilities, he de¬ 
termined scaling laws for the numerical viscosity and resistivity of 
the code as a function of resolution and initial conditions. Given 
our numerical setup, the numerical viscosity and resistivity of the 
code are <10^ cm^s"^ if the relevant length scale is covered by at 
least 20 zones (see also section 4.2.5.1). 

To investigate the nature of the parasitic instabilities that are 
responsible for quenching the MRI, we used both Fourier analysis 
and the local magnetic held structure during MRI termination. In 
models #2 to #11, we And parasitic instabilities of KH type devel¬ 
oping at an angle consistent with 0p = 45° independent of resolu¬ 
tion. The horizontal wavelength of the parasitic modes ranges from 
Tp ^ 0.50 km in the coarsest resolved model #2 to dp ^ 0.23 km in 
the best resolved models #10 and #11. 

Remarkably KH instabilities develop even in model #2, al¬ 
though it was simulated with the coarsest grid which noticeably 
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affects the structure of the channel modes. Their identification was 
not straightforward, however. From the patterns recognizable in the 
magnetic field structure during MRI termination, one can exclude 
TM and other obvious numerical artefacts, e.g., arising from im¬ 
perfect boundary conditions. We also find patterns resembling un¬ 
derresolved vortex rolls which we attribute to KH modes, although 
without the help of Fourier analysis, our experience gained from 
observing similar patterns in better resolved simulations, and most 
importantly a theoretical model (GX94; [Pessah|2010] ), this identifi¬ 
cation would have been impossible. 

The value of the Maxwell stress at MRI termination de¬ 
pends on the grid resolution (see Tab.[^ and green circle symbols 
in Fig.[^. It is larger for coarse grids, the only exception be¬ 
ing model #2. Using the MP9 reconstruction scheme, the Maxwell 
stress converges to a value of Mrtp - 0.73 x 10^^ in models 
#10 and #11 simulated with 67 and 134 grid cells per MRI chan¬ 
nel, respectively. The dashed line in Fig. 13 shows that Mrcp follows 
roughly a power law with index -1/3 as a function of the number 
of zones per MRI channel, except for the under-resolved model #2 
with only eight zones per channel, in which MRI growth terminates 
at a much smaller value of Mrcp. 

This behaviour of Mrcp can be explained within the para¬ 
site model. The KH modes responsible for MRI termination grow 
from the shear layers created by the MRI channels. During their 
growth, these secondary instabilities develop structures that are 
even smaller than the MRI channels. Therefore, to resolve these 
structures a finer grid is necessary than for the MRI channels them¬ 
selves. As a result, simulations may not be converged even if the 
numerically obtained MRI growth rate is close to the theoretical 
value. Because badly resolved KH modes grow at a rate below the 
theoretical expectation, parasitic instabilities reach the necessary 
amplitude to disrupt the channels slightly later during the evolu¬ 
tion, and the MRI will terminate at a somewhat higher amplitude 
than in well-resolved simulations. 

The appearance of the outlier (model #2) in Fig.[^is not sur¬ 
prising. We have already seen that MRI channels are underresolved 
in this model. Therefore, one would not expect to find any rea¬ 
sonable termination amplitude for this model, i.e. it is probably a 
chance coincidence that this model gives roughly the same results 
as the other better resolved models. 

For completeness, we included also models #12—#18 in Fig.p~3] 
(red diamond symbols), for which the parasitic modes are restricted 
not only by the grid resolution but also by the box size and some 
physical viscosity. The data of these models lie all above the dashed 
line indicating the infiuence of resolution alone. Finally, the asterisk 
symbol marks the result of the ideal MHD model #1, which we 
discuss in the next subsection. 


4.2.5 Influence of viscosity and resistivity 

4.2.5.1 Viscosity We have performed the reference simulation 
(model #7; see Tab.[T]), and the set of models described above with¬ 
out any explicit physical viscosity. To study whether a finite physi¬ 
cal viscosity corresponding to Re ~ 100 (for which the equations of 
ideal hydrodynamics should still approximately hold) leads to more 
than mere quantitative changes, we compare the results of models 
#9 and #12. As expected, in both of these two models, which differ 
only in the value of the hydrodynamic Reynolds number (Re = ^ 
for model #9, and Re ~ 100 for model #12) the growth of the MRI 
is terminated by KH modes. 

Even though no physical viscosity was used for model #9, it 
is affected by a non-zero numerical viscosity. To infer the quantita- 
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Figure 13. MRI termination amplitudes as a function of grid resolution 
for models #1 (blue asterisk), #2^11 (green circles), and #12^18 (red 
diamonds), respectively. The dashed line represents a power law in the grid 
resolution with an index of -1/3. 
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five importance of the numerical viscosity, we first compare models 
which differ only in grid resolution in the inviscid limit, e.g., model 
#9, where the evolution is still notably affected by numerical vis¬ 
cosity, and models #10 and #11, which are already converged, as 
we have argued in the previous section. We note that the numer¬ 
ical viscosity accounts for the ^ 22% difference in the volume- 
averaged Maxwell stress at termination (i.e. the parasitic instabili¬ 
ties are somewhat under-resolved), and has a negligible impact on 
the MRI growth rate. 

We now turn to the comparison of an inviscid model (#9) with 
one where a physical viscosity has been included (#12). Quanti¬ 
fying the exact contribution of the numerical viscosity in a model 
that includes also physical viscosity is much more involved. The 
reasons for this difficulty are multiple. On the one hand, the numer¬ 
ical viscosity is a scale dependent effect, affecting more the smallest 
scales. Hence, it is relatively more important for the development 
of parasitic modes than for the MRI growth. On the other hand, 
the volume-averaged Maxwell stress at termination depends chiefly 
on a delicate balance between two competing factors. First, models 
computed with a finite physical viscosity have a smaller yMRi- Thus, 
for a similar growth time, Mrcp should be smaller than for models 
computed without a physical viscosity. Secondly, this effect is off¬ 
set, however, by the also smaller growth rate of the KH parasitic 
modes, which defers the saturation process and allows for a further 
growth of Mrcf). Thereby, our comparison is of qualitative nature 
only. In view of the former reasoning, we tentatively attribute the 
~ 4% difference in A\r(p between the inviscid model #9 and the vis¬ 
cous model #12 to the effect of a finite viscosity in the latter one, 
though we cannot unambiguously say which part of this difference 
is caused by numerical viscosity. To support the argument that a 
finite physical viscosity and not a numerical viscosity is driving the 
dynamics, we note that the MRI growth rate decreases significantly 








































18 T. Rembiasz et al. 


in model #12 (as also in all other viscous models from #13 to #18). 
This expected behaviour (see Section \22\ implies that including 
a physical viscosity causes measurable differences already in the 
MRI growth phase. 

4.2.5.2 Resistivity In a final set of simulations, we investigated 
the influence of a finite (physical) resistivity on MRI termination. 
Model #1, simulated in ideal MHD, represents the typical condi¬ 
tions prevailing in core collapse supernovae close to the surface of 
PNS, whereas in models #19 and #20 we explored the (opposite) 
limit of a very high resistivity. With the latter simulations we also 
address a prediction of|P essah| ( [20T0] ), namely that for Rm < T the 
MRI should be terminated by TM. We performed the ideal MHD 
simulation of model#! in the default size box, the results being 
very similar to those obtained for model#?, which we simulated 
with the identical setup, but with ^ 100. 

As expected, the MRI grew in the ideal MHD simulations at a 
somewhat higher rate (yMRi = 1137 s"^, because its growth was not 
damped by resistivity. We found that the MRI is terminated by par¬ 
asitic KH instabilities in this model, too. The parasites develop at 
an angle 0p ^ 45° in accordance with the parasite model of GX94. 
The Maxwell stress at MRI termination is Mrcp - 0.96 x 10^^ G^, a 
value which is similar to those obtained for models #7 and #12 (the 
latter one being simulated with ^ = 100). From these find¬ 

ings we conclude that for Reynolds numbers Rq^R^ ^ 100 viscos¬ 
ity and resistivity have no strong influence on the process of MRI 
termination and the saturation value of the Maxwell stress. How¬ 
ever, one should not forget that in all three simulations (i.e. mod¬ 
els #1, #19, and #20), numerical resistivity and numerical viscosity 
do contribute to the ’total’ hydrodynamic and magnetic Reynolds 
numbers (by lowering them). 

In the limit of very high resistivity, i.e. for R^ 1, the MRI 
modes that grow fastest in ideal MHD are completely suppressed 
by magnetic dissipation (PC08). Only modes with sufficiently long 
wavelengths survive the interplay between MRI amplification and 
resistive damping in this limiting case, because magnetic dissi¬ 
pation is weaker for these modes. The wavelength of the fastest- 
growing mode is Tmri oc (PC08). To keep Tmri ~ 0.333 km 

in models #19 {R^ - 0.1) and #20 {R^ - 0.05), we had to reduce 
the initial magnetic field strength to values of = 3.25 x 10^^ G 
and boz - 1.63 x 10^^ G , respectively. 

We find that the fastest MRI mode grows at a considerably 
slower rate than in the ideal MHD limit in both models, the mea¬ 
sured growth rates being yMRi = 80 s"^ (model #19) and yMRi = 
34 s“^ (model #20). These values agree very well with those of 
PC08, who predicted 7 mri = 79 s~^ and yMRi = 40 s"^ respectively. 
PC08 also predicted that for Re ^ 1 and R^ - 0.1, we should find 
0v = 2° and 0^ = 94°, i.e. the velocity channels and the magnetic 
field channels should almost be aligned with the r axis and the 0 
axis, respectively. Indeed, our simulations show that \b(f,\ » \brl 
i.e. 0^ ^ 90°, during the exponential growth phase (see Fig. [E}. 

Finally, according to |Pessah| ( |2010| ), the MRI should be termi¬ 
nated by TM and not KH instabilities for < 1- For R^ = 0.1, the 
dominant parasitic mode should have a wavelength 

dp = Ttm — 2 .1 /Imri (61) 

developing at an angle 0p = 0^? = 94°. A Fourier analysis of the 
MRI modes and their parasites (see Fig.[^ shows that shortly be¬ 
fore MRI termination (t ^ 160 ms) the energy stored in horizontal 
Fourier modes increases significantly growing super-exponentially 
with time. However, because of numerical noise introduced by our 
imperfect boundary conditions, we could neither identify parasitic 



time [ms] 

Figure 14. Same as Fig.|ll|but for model #19 simulated with a high resistiv¬ 
ity (Rm = 0.1). Note the very different time-scale at which the MRI evolves 
in this case. 

instabilities nor determine their angle with the help of the Fourier 
analysis. Therefore, the type of the parasitic instability terminating 
the MRI in model #19 (as well as in model #20) remains unknown. 
It is probably neither KH nor TM, because we could not recognise 
the characteristic features of none of these parasites. The bottom 
line is that based on models #19 and #20, we cannot confirm the 
predictions of |Pessah| ( |2010| ) that for R^ < 1 the MRI is termi¬ 
nated by TM. As the highly resistive limit is of no direct relevance 
for core collapse supernovae, we did not investigate this regime in 
more detail. 


5 SUMMARY AND DISCUSSION 

[Akiyama et al.| ( |2003| ) first pointed out that the MRI can develop 
in CCSNe and could lead to a strong amplification of the mag¬ 
netic field during the proto-neutron star phase. The presence of a 
dynamically relevant magnetic field is of great interest because it 
may have implications for the supernova explosion mechanism and 
for explaining the origin of neutron stars magnetic fields. However, 
for the typical magnetic field strength present in supernova pro¬ 
genitors, it is necessary to resolve numerically MRI channels with 
a length scale ~ 1 m, which makes global simulations unfeasible 
even with present day supercomputers. To overcome the limita¬ 
tion of global simulations |Obergaulinger et ah] ( |2009| ) performed 
semi-global simulations in a kilometre size box located near the 
surface of the proto-neutron star. They confirmed the hypothesis of 
[Akiyama et al.| ( [200^ showing that the MRI can strongly amplify 
the magnetic field and parasitic instabilities are able to destroy the 
MRI channels, thereby quenching the amplification of the magnetic 
field. However, they were unable to identify the agent terminating 
the MRI growth. Hence, [Obergaulinger et al.| ( [2009| ) could not give 
an upper limit on the MRI-driven magnetic field amplification. 
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This work tries to shed some light on the termination process. 
We numerically studied the evolution of the MRI in the context 
of a CCSN. To this end, we performed a set of local 3D resistive- 
viscous MHD simulations in a box representing typical conditions 
in a proto-neutron star. In all of our simulations we observed the 
growth of the MRI and its subsequent termination. We identified 
that secondary KH instabilities are responsible for MRI termination 
by acting as parasites on the MRI channels. Our results are consis¬ 
tent with the predictions of the parasite model proposed by GX94 
and further developed by [Latter et al.| ( [200^ and |Pessah| ( |2010| ). 
Hence, the parasite model, which is based on a local stability anal¬ 
ysis of the system, allows some insight into the termination process 
and provides a valuable guidance for interpreting numerical results. 
This conclusion was obtained by an analysis of our numerical sim¬ 
ulations based on MRI theory and the parasite model in a number 
of tests: 

(i) We observe an exponential growth of MRI channels of a size 

compatible with the predicted value, Tmri, for the fastest-growing 
mode (BH91). The growth rates measured in the numerical simu¬ 
lations agree with those obtained from the local analysis within the 
simulation box. The presence of viscosity or resistivity lowers the 
values of the growth rates as expected from theory. In the limit of 
large hydrodynamic and magnetic Reynolds numbers, ^ L 

the angle between velocity (magnetic field) in the MRI channels 
and the radial direction is close to the theoretically expected one in 
the ideal MHD limit by GX94, i.e. 0^ = 45° (0^, = 135°). 

(ii) For simulations with large hydrodynamic and magnetic 
Reynolds numbers, R^.R^ ^ L the MRI is terminated by KH in¬ 
stabilities, as predicted by GX94 and |Pessah| ( |2010| ) for this regime. 
To identify the KH instabilities we had to reduce the time interval 
between outputs of our simulations around termination to a value 
below the typical growth time of the parasites, in our case 0.1 ms. 
This explains why jObergaulinger et al.| ( [200^ did not observe the 
development of KH instabilities, since they employed insufficiently 
frequent output. To help the analysis we project the magnetic field 
components into the directions 0v and 0^. Shortly before termina¬ 
tion we observe the development of vortex rolls at the position of 
the shear layer located along an angle 0v = 45°. We identify the vor¬ 
tices as a consequence of KH instabilities, which eventually disrupt 
the channels causing turbulence. This behaviour is in agreement 
with the predictions of GX94 and |Pessah| < |2010| ). 

(iii) Analysing the numerical simulations by means of Fourier 
transforms, we have been able to determine the properties of the 
parasitic instabilities, i.e. their horizontal wavelengths and the an¬ 
gles at which they develop. Our results confirm the theoretical pre¬ 
dictions of GX94 and |Pessah| ( |2010| ) that for R^^R^ ^ 1 the domi¬ 
nant parasitic instabilities are KH modes developing parallel to the 
MRI velocity channels (i.e. at an angle 0p ^ 45° in the (r, 0) plane). 
The horizontal length of these KH modes is in a reasonable agree¬ 
ment (within a factor of 2) with the predictions of|P ess SIl poTot . 

(iv) Motivated by the good agreement with the parasite model, 
we explored the regime of very high resistivity (i.e. R^ < 1), al¬ 
though it is not of direct relevance for CCSNe. In this regime [Pes^ 
|sah| ( |20I0] ) predicts that the MRI should be terminated by TMs. To 
this end we performed two simulations with = 0.1 and 0.05, but 
we were unable to identify the agent terminating the MRI growth, 
because no clear signature of TM or KH instability could be iden¬ 
tified in our simulations. We conclude that possibly even higher re¬ 
sistivities are required for TM to become the dominant secondary 
instability. A possible explanation of this discrepancy between the¬ 
ory and our simulations is the fact that the calculations of jPessah] 


( |2010| ) are based on somewhat simplified assumptions. Therefore, 
his results should be treated more like a guideline rather than exact 
predictions. 

To confirm the above conclusions, we have studied systemati¬ 
cally the effect of resistivity, viscosity, box size, and grid resolution. 
There are a number of numerical artefacts that can affect the simu¬ 
lations: 

(i) In axisymmetric (2D) simulations, the MRI growth is termi¬ 
nated by TM. This was already reported by jObergaulinger et ah] 
( |2009| ), who observed TM in their 2D ideal MHD simulations, 
in accordance with theoretical results obtained by |Pessah| ( |2010| ). 
jObergaulinger et aL] ( |2009j l argued that this resistive MHD instabil¬ 
ity must have been triggered by numerical resistivity. In 2D simula¬ 
tions only axisymmetric parasitic modes can develop. As axisym¬ 
metric KH modes are strongly suppressed by the magnetic field 
tension of the MRI channels, TM become the dominant secondary 
instability. They suffer less strongly from this constraint, develop 
faster than axisymmetric KH modes (but slower than KH modes 
would grow in full 3D), and terminate the MRI growth. As a re¬ 
sult the MRI in 2D is terminated at unrealistically large magnetic 
stresses and continues to grow after termination, a behaviour not 
observed in our 3D simulations. 

(ii) We performed a set of simulations to study the dependence 
of the properties of the parasite modes. Using a 9th-order spatial 
reconstruction scheme (MP9), we reached convergent results only 
in simulations with at least 60 zones per MRI channel, i.e. the value 
of the magnetic stress at termination differed by less than 10% be¬ 
tween the two highest resolution runs. This grid resolution is signif¬ 
icantly higher than the one necessary to obtain convergence in the 
growth rate of the MRI, for which 8 zones are sufficient to obtain 
a 10% accuracy. We note that the required number of zones will 
be significantly higher, if lower order reconstruction methods are 
used. Our result is not surprising when viewed in the light of the 
parasite model. The KH instability is triggered by the shear layer 
between MRI channels. At this layer structures develop that are 
much finer than the width of the channel itself. Failing to capture 
the KH instabilities properly because of a lack of resolution leads 
to artificially large magnetic stresses at MRI termination. However, 
the qualitative behaviour of the flow seems not to be affected by 
a lack of resolution. Even for our lowest resolution simulation (8 
zones per channel) we are able to identify KH instabilities develop¬ 
ing at 0p ^45° as the main termination agent. 

(iii) We studied how the box size can affect the development 
of parasitic instabilities. In 3D simulations with an azimuthal box 
length of at least ^ 1/1 mri the MRI is terminated by KH modes, 
whereas simulations with an azimuthal box length < ITmri gave 
very similar results as 2D simulations and the MRI was terminated 
by TM (because KH modes are suppressed; compare models #16 
and #17 in Tab.[^. Taking into account that the parasitic KH insta¬ 
bilities develop at a 45° angle, one should use a box size of at least 
1 ~ 4 mri is in the radial direction. However, we recommend us¬ 
ing larger boxes in the horizontal directions to reduce the influence 
of the boundary conditions on the development of the parasitic in¬ 
stabilities. In the vertical direction it is sufficient to use a box size 
of Tmri, i.e. it is sufficient to consider the evolution of one single 
MRI channel to capture the termination process correctly. Deter¬ 
mining the minimum box size has been critical to be able to per¬ 
form the simulations with the highest resolution presented in this 
work (model #11 with 134 zones per channel). 

Our results have some important implications for the commu- 
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nity modelling magnetorotational collapse of stellar cores. Study¬ 
ing the termination of the exponential growth phase of the MRI, we 
learned that it involves parasitic KH instabilities, but does neither 
depend on the physical resistivity nor viscosity present in CCSNe. 
The effects of the interaction of neutrinos and matter, which at dif¬ 
ferent locations in the PNS can be described either as an effective 
viscosity or by a drag term, affect the properties of the MRI, in 
particular its growth rate and wavelength ( [Guilet et al.|2015| ). Con¬ 
sequently, this interaction should induce modifications of the termi¬ 
nation amplitude. On the other hand, we do not expect these effects 
to change the most important qualitative result of our study, namely 
that the MRI is not terminated by TM, but by the KH instability. 

The presence of background flows in CCSNe (like convection 
or turbulence) can affect the growth and termination of the MRI, 
too. In 2D global simulations, |Cerda-Duran et aL] ( |2008| ) observed 
the growth of coherent channel flows, albeit deformed, while|Saw^ 
|et al.| ( [20T^ observed the disruption of channel flows on timescales 
of milliseconds due to the presence of a dynamical background. 
Both, the parasitic instabilities and the effect of the background 
flow, seem to disrupt channel flows on similar timescales. There¬ 
fore, it is likely that both effects will be of relevance in studying the 
MRI termination in global simulations. 

Furthermore, the presence of parasitic instabilities implies the 
existence of a maximum magnetic stress at termination. This limits 
the ability of the MRI to amplify the magnetic held to dynami¬ 
cally relevant values in CCSNe. The turbulence triggered by the 
MRI may amplify the global magnetic held further, if conditions 
for dynamo action are encountered. A study of the latter process 
is beyond the scope of this work, and will be addressed in a future 
publication. The MRI has been traditionally invoked in the CCSN 
community as a justification to use supernova progenitors with an 
artificially enhanced initial magnetic held strength, 10^-10^ times 
larger than the values expected in stellar evolution models. This ar¬ 
gumentation should be revisited if we want to move towards a more 
realistic modeling of magnetorotational core collapse. 

Unfortunately, given the minimum numerical requirements 
that we provide in this work, 3D simulations seem nowadays 
unfeasible. To resolve MRI channels of 1 m length scale in a 
proto-neutron star of about 30 km radius, with a resolution of 
~ 60 zones per channel, it would require about 5 x 10^^ grid 
zones. Evolving this system for a typical dynamical timescale of 
the proto-neutron star evolution, ~ 100 ms, would require about 
10 ^ time steps, several zettabytes of memory and would con¬ 
sume ~ 10^^ CPU hours (10^ iterations times 5 x 10^^ zones times 
100 operations/zone at SOOGflops). With current petascale super¬ 
computers this would amount to ~ 3 X 10^ yr of uninterrupted com¬ 
putation with 100 000 cores. A simulation which includes the whole 
iron core with a radius of ~ 1000 km would be even more demand¬ 
ing. Reducing the dimensionality of the system by imposing ax- 
isymmetry reduces the computational time by a factor ~ 4 x 10^. 
The use of adaptive mesh refinement (AMR) techniques would al¬ 
leviate, but not solve the problem. 

The CCSN community has been performing 2D and 3D simu¬ 
lations with resolutions of up to 12.5 m (2D; |Sawai et al.|20T^ and 
15.6 m (3D; [Masada et al.|2015| l. We note that these authors could 
afford such high resolutions only by accepting other limitations in 
terms of geometry and input physics, while other simulations use 


resolutions of several 100m (Obergaulinger et al.|2006a|b[ 

Shibata 

et al.|20061|Burrows et al.|20071|Cerda-Duran et al.|2008 

jSchei- 


degger et al.||2008| I Mikami et al.||2008[ IKuroda & Umeda|2010 

Takiwaki & KotakepOll jWinteler et al.|2012[|Mosta et al.|2014 


Obergaulinger et al.|2014 1, which even for the artificially enhanced 


initial magnetic held, fail to resolve the termination process of the 
MRI. Given the above computational requirements, the community 
should rethink the way to model this scenario. The agreement that 
we have found here between local simulations and the parasitic in¬ 
stability analysis of GX94 and jPessahj ( |2010j ) make us hope that 
there could be a way out. If we can understand how MRI and tur¬ 
bulence work at sub-meter scales, this information could be incor¬ 
porated in global simulations using appropriate subgrid models. As 
a first step along this path we plan to gain a deeper understanding 
in the maximum magnetic stress achievable by the MRI and in the 
properties of its turbulent saturated state. 
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